#include <iostream>
#include <complex>
#include "bessel.h"

using namespace std;

int main()
{
    cout << "Bessel function check.\n";

    complex <double> z, nu, ans;

    int KODE, NZ, IERR;

    KODE = 1;

    nu = complex <double> (1000.0, -50.0);
    z= complex<double> (1000.0, 0.0);


    cout << nu << " \t" << z  << "\n"
	 << "J : " << BesselJ(nu, z, KODE, &NZ, &IERR) << "\t"
	 << NZ << "\t" << IERR << "\n"
	 <<  "Y : " << BesselY(nu, z, KODE, &NZ, &IERR) << "\t"
	 << NZ << "\t" << IERR << "\n"
         << "H1 : " << HankelH1(nu, z, KODE, &NZ, &IERR) << "\t"
	 << NZ << "\t" << IERR << "\n"
	 << "H2 : " << HankelH2(nu, z, KODE, &NZ, &IERR) << "\t"
	 << NZ << "\t" << IERR << "\n"
	 << "DJ : " << DBesselJ(nu, z, KODE, &NZ, &IERR) << "\t"
	 << NZ << "\t" << IERR << "\n"
	 <<  "DY : " << DBesselY(nu, z, KODE, &NZ, &IERR) << "\t"
	 << NZ << "\t" << IERR << "\n"
         << "DH1 : " << DHankelH1(nu, z, KODE, &NZ, &IERR) << "\t"
	 << NZ << "\t" << IERR << "\n"
	 << "DH2 : " << DHankelH2(nu, z, KODE, &NZ, &IERR) << "\t"
	 << NZ << "\t" << IERR << "\n";

	 ;



    return 0;
}

/* Benchmark results

Bessel function check.
(1000,-50) 	(1000,0)
J : (-0.458395,38.1709)	0	0
Y : (-38.1709,-0.458376)	0	0
H1 : (-1.86448e-05,-1.86397e-05)	0	0
H2 : (-0.916771,76.3418)	0	0
DJ : (8.34551,8.53214)	0	0
DY : (-8.53214,8.3455)	0	0
DH1 : (8.43297e-06,-1.91299e-07)	0	0
DH2 : (16.691,17.0643)	0	0

*/
